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Turbulent suspensions of heavy particles in incompressible flows have gained much attention in 
recent years. A large amount of work focused on the impact that the inertia and the dissipative 
dynamics of the particles have on their dynamical and statistical properties. Substantial progress 
followed from the study of suspensions in model flows which, although much simpler, reproduce 
most of the important mechanisms observed in real turbulence. This paper presents recent de- 
velopments made on the relative motion of a pair of particles suspended in time-uncorrelated and 
spatially self-similar Gaussian flows. This review is complemented by new results. By introducing a 
time-dependent Stokes number, it is demonstrated that inertial particle relative dispersion recovers 
asymptotically Richardson's diffusion associated to simple tracers. A perturbative (homogeneiza- 
tion) technique is used in the small-Stokes-number asymptotics and leads to interpreting flrst-order 
corrections to tracer dynamics in terms of an effective drift. This expansion implies that the corre- 
lation dimension deflcit behaves linearly as a function of the Stokes number. The validity and the 
accuracy of this prediction is confirmed by numerical simulations. 
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I. INTRODUCTION 

The current understanding of passive turbulent trans- 
port profited significantly from studies of the advection 
by random fields. In particular, fiows belonging to the 
so-called Kraichnan ensemble — i. e. spatially self-similar 
Gaussian velocity fields with no time correlation — which 
was first introduced in the late 1960's by R.H. Kraich- 
nan led in the mid 1990's to a first analytical de- 
scription of anomalous scaling in turbulence (see for a 
review). More recently, much work is devoted to a gener- 
alization of this passive advection to heavy particles that, 
conversely to tracers, do not follow the flow exactly but 
lag behind it due to their inertia. The particle dynamics 
is thus dissipative even if the carrier flow is incompress- 
ible. This paper provides an overview of several recent 
results on the dynamics of very heavy particles suspended 
in random flows belonging to the Kraichnan ensemble. 

The recent shift of focus to the transport of heavy par- 
ticles is motivated by the fact that in many natural and 
industrial flows flnite-size and mass effects of the sus- 
pended particles cannot be neglected. Important applica- 
tions encompass rain formation [1, H 01 and suspensions 
of biological organisms in the ocean [6|,l2l, Q . For practical 
purposes, the formation of particle clusters due to inertia 
is of central importance as the presence of such inhomo- 
geneities significantly enhances interactions between the 



suspended particles. However, detailed and reliable pre- 
dictions on collision or reaction rates, which arc crucial 
to many applications, are still missing. 

Two mechanisms compete in the formation of clusters. 
First, particles much denser than the fiuid are ejected 
from the eddies of the carrier fiow and concentrate in the 
strain-dominated regions 0. Second, the dissipative dy- 
namics leads the particle trajectories to converge onto a 
fractal, dynamically evolving attractor p^.[ll|. In many 
studies, a carrier velocity field with no time correlation — 
and thus no persistent structures — is used to isolate the 
latter effect. As interactions between three or more par- 
ticles are usually sub-dominant, most of the interesting 
features of mono-disperse suspensions can be captured by 
focusing on the relative motion of two particles separated 
by R. 
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where dots denote time derivatives and r the particle re- 
sponse time. The fluid velocity difference Su is a Gaus- 
sian vector field with correlation 

{Su\r,t)Su^{r',t')) = 2h'^ {r - r') 5{t - t'). (2) 

In order to model turbulent fiows, the tensorial struc- 
ture of the spatial correlation b^^ (r) is chosen to ensure 
incompressibility, isotropy and scale invariance, namely 

b'^r) = Di r^''[{d -l + 2h) S'^ - 2hr'r^ /r% (3) 

where h relates to the Holder exponent of the fiuid veloc- 
ity field and Di measures the intensity of its fluctuations. 
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In particular, h = 1 corresponds to a spatially difFeren- 
tiable velocity field, mimicking the dissipative range of 
a turbulent flow, while h < 1 models rough flows, as in 
the inertial range of turbulence. In this paper we mostly 
focus on space dimensions d = 1 and d — 2; extensions 
to higher dimensions are just sketched. 

The above depicted model flow has the advantage that 
the particle dynamics is a Markov process. In particular, 
Gaussianity and ^-correlation in time of the fluid velocity 
field imply that the probability density p{r, v, t\ro, Vq, to) 
of finding the particles at separation R{t) = r and with 
relative velocity R{t) — v a.t time t, when R{to) — Tq and 
R{to) = is a solution of the Fokker-Planck equation 

dtp + T.(^r- Idl) («» - E ^ ^vdiP = 0, (4) 

with the initial condition p{r,v, to) =^ S{r — ro) d{v — Vq). 
To maintain a statistical steady state, the Fokker-Planck 
equation ^ as well as the stochastic differential equa- 
tion (III) should be supplemented by boundary conditions, 
here chosen to be reflective at a given distance L. 

For smooth flows {h = 1), the intensity of inertia is gen- 
erally measured by the Stokes number St, defined as the 
ratio between the particle response time t and the fluid 
characteristic time scale. For St ^ 0, particles recover the 
incompressible dynamics of tracers. In the opposite limit 
where St is very large, inertia effects dominate and the 
dynamics approaches that of free particles. In the above 
depicted model, the Stokes number is defined by non- 
dimensionalizing r by the typical fluid velocity gradient, 
i.e. St = DiT. Note that by rescaling the physical time 
by T, it is straightforward to recognize that the dynamics 
depends solely on St. 

Similarly it can be checked that in rough flows {h<l) 
— with an additional rescaling of the distances by a fac- 
tor (_Dit)^/(^^^''^ — the dynamics of a particle pair at 
a distance r only depends on the local Stokes number 
St(r) — DiT/r^'^^~'^\ This dimcnsionless quantity, first 
introduced in 12] and later used in [l^, is a general- 
ization of the Stokes number to cases in which the fluid 
turnover times depend on the observation scale. At large 
scales, St(r) — > and inertia becomes negligible. Parti- 
cle dynamics thus approaches that of tracers. At small 
scales, St(r) — > oo and the particle and fluid motions 
decorrelate, so that the inertial particles move ballisti- 
cally. In both the large and small Stokes number asymp- 
totics, particles distribute uniformly in space, while inho- 
mogeneities are expected at intermediate values of St(r). 

The paper is organized as follows. In Section|TTl an ap- 
proach originally proposed in is used to reduce the 
dynamics of the particle separation to a system of three 
stochastic equations with additive noises. This formula- 
tion is useful for both numerical and analytical purposes, 
particularly when studying the statistical properties of 
particle pairs. In Section lllll we introduce the correlation 
dimension to quantify clustering as well as the approach- 
ing rate which measures collisions. Numerical results for 



these quantities are reported. In Sect ion [TVl we introduce 
the notion of time-dependent Stokes number which makes 
particularly transparent the interpretation of the behav- 
ior of the long-time separation between particles. We 
show how Richardson dispersion, as for tracers, is recov- 
ered in the long time asymptotics. Section |V] briefly sum- 
marizes some exact results that can be obtained for the 
one-dimensional case. Sections I VII and fVIII are dedicated 
to the small and large Stokes number asymptotics, re- 
spectively. In particular, the former one presents an orig- 
inal perturbative approach which turned out to predict, 
in agreement with numerical computations, the behavior 
of the correlation dimension that characterizes particle 
clusters. Finally, Section rVIIII encompasses conclusions, 
open questions and discusses the relevance of the consid- 
ered model for real suspensions in turbulent flows. 

II. REDUCED DYNAMICS FOR THE 
TWO-POINT MOTION 

In this Section we focus on planar suspensions {d = 2). 
Following the approach proposed in [l4| and with the 
notation R ^ \R\, the change of variables 

ai = {L/R)^+''R- R/L^, (5) 
a2 = {L/R)^+''\R/\R\/L\ (6) 
P = (R/L)'-'^ (7) 

is introduced to reduce the original system oi 2d = A 
stochastic equations to the following one of only three 
equations 

ai = -(Ji/T-[hal-aj]/p + VCr]i, (8) 

(T2 = -Ta/r - (1 + h)aia2/p + + 2h)C r]2, (9) 
p=il-h)ai, (10) 

where C — 2Di/ {tL^~^)^ and rji denote two independent 
standard white noises. Reflective boundary conditions 
at i? = L in physical space imply reflection at p = 1. 
Note that txi and (T2 are proportional to the longitudinal 
and to the transversal relative velocities between the two 
particles. In the smooth case {h = 1), one has p = 1 and 
equations ([8]) and ([9]) decouple from pO)l . The particle 
separation R then evolves as 

R = ai{t)R. (11) 

Besides this simple evolution and the reduction of the 
number of variables from 2d to only three, the change 
of variables {R,R\ ^ {p, 0-1,(72} has several other ad- 
vantages. For instance the noise, which is multiplica- 
tive in the original dynamics ([T]), becomes additive in 
the reduced system ([5t - ([TU)) . However, this simplifica- 
tion is counter-balanced by the presence of nonlinear drift 
terms. Note that in dimensions higher than two, there 
is an additional term cx l/cr2, which is due to the Ito 
formula [ll,[3. 
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FIG. 1: Sketch of the reduced dynamics ((U - ^ for h = 0.7. 
The dotted lines represent the drift. The sohd line depicts a 
random trajectory with St(L) = 1. (a) full (ai, a"2, p)-space, 
(b) projection on p = plane, and (c) on the a2 = plane. 

Figure[l] sketches the deterministic drift and shows a 
typical trajectory in the reduced space. This dynam- 
ics can be qualitatively described as follows. The line 
CTi =(T2 = acts as a stable fixed line for the drift. Hence 
a typical trajectory spends a long time diffusing around 
it, until the noise realization becomes strong enough to 
let the trajectory escape from the vicinity of this line. 
Whenever this happens with a positive longitudinal rel- 
ative velocity (cri > 0), the trajectory is pulled back to 
the stable line by the quadratic terms in the drift. Con- 
versely, if CTi < and /icr^ +(71^— (t| < 0, the drift pushes 
the trajectory towards larger negative values of ai. Then 
the particles get closer to each other and p decreases, 
until the quadratic terms in equations ^ and ^ be- 
come dominant. The trajectory then loops back in the 
(<Ti, cr2)-plane, approaching the stable line from its right. 
It is during these loops that the inter-particle distance R 
becomes substantially small. The loops thus provide the 
main mechanisms for cluster formation. 



Velocity statistics 

Numerical simulations show that the probability den- 
sity function (pdf ) of the longitudinal relative velocity cri 
displays algebraic tails at large positive and negative val- 
ues (see Fig. [5]). As will become clear in the sequel, these 
power-law tails are a signature of the above-mentioned 
large loops. Let us consider the cumulative probability 
P^{(j) = Pr (fTi < a) for cr<C— 1. This quantity can be 
estimated as the product of (i) the probability to start a 
sufficiently large loop in the (cri, cr2)-plane that reaches 
values smaller than a and (ii) the fraction of time spent 
by the trajectory at ai < a. Within a distance of the 
order of unity from the line cri = (T2 = 0, the quadratic 
terms in the drift are subdominant and can be disre- 
garded. Then cri and ai can be approximated by two 



FIG. 2: Log-log plot of the pdf of en for St(L) = 1 for five 
values of the fluid Holder exponent h. Power-law tails are 
always observed, p(cr) cx |cr|~". Inset: exponent a versus ft; 
the dashed line is the theoretical prediction a = 1-1-2/ ft. 

independent Ornstein-Uhlenbeck processes. Conversely, 
at sufficiently large distances from that line, only the 
quadratic terms in the drift contribute and the noise is 
negligible. 

Within this simplified dynamics, a loop is initiated at a 
time to for which a\(t^<—\ and cr2 (to ) <C | ci (t o ) | . Once 
these conditions are fulfilled, the trajectory performs a 
loop in the (cri, cr2)-plane and both |cri(t)| and cr2(<) be- 
come very large. The maximum distance from the stable 
line, which gives an estimate of the loop radius, is reached 
when cr2 is of the order of \(J\\. Let t* denote the time 
when this happens, i.e. cr2(t*)/|cri(t*)| = 0(1). When 
neglecting the noise, this condition leads to the following 
estimate for the loop radius 

|ai(r)| (X [ai(to)+p(io)/r] ^^1(^0)1" (Ta2(to))-' , (12) 

see [l3| for details. In order to reach velocity differences 
such that cTi < cr<C — 1, the radius of the loop has to be 
larger than |cr|. From this implies that cr2(to) has 
to be smaller than |cr|~^/''. In order to evaluate contri- 
bution (i), one has to estimate the probability to have 
o'i(to) ^ ^1 and 0-2(^0) < |cr|~^/'' from the dynamics 
in the vicinity of the origin. Approximating the two ve- 
locity differences cri and 02 by independent Ornstein- 
Uhlenbeck processes close to the line ai = cr2 = 0, the 
first condition gives an order-unity contribution, while 
the second has a probability c>c jcrj"-'^/''. For estimating 
(ii), we neglect the noise in the dynamics far from the 
stable line. The probability is then given by the frac- 
tion of time spent at cri < cr which is proportional to 
cr2(to) cx |cr|^^/''. Put together, the two contributions 
yield P'^{x) oc |cr|^^/'' when cr <C —1. Thus the negative 
tail of the pdf of cri behaves as oc |cr|~", with a = l-|-2/ft. 

During the large loops, the trajectories equally reach 
large positive values of cri and of cr2. Again the fraction 
of time spent at both cri and cr2 larger than cr 3> 1 can be 
estimated as . Hence, the pdf of both longitudinal 

cri and transversal cr2 velocity differences have algebraic 
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left and right tails with exponent a. Both tails are de- 
picted in Fig. [21 where the inset shows that the numerical 
measurements are in good agreement with the predicted 
value of a. The relation between a and the Holder ex- 
ponent h implies in particular that a = 3 in the smooth 
case, while it increases with decreasing h. Moreover, it 
follows straightforwardly from ([5]) - PH)) that during the 
loops p{t) IX. p{to)'^ when p{to) ^ 1. Hence it becomes 
less and less probable to reach smaller values of p as ft. 
decreases. In other words, particle clustering should be 
very strong for smooth flows and becomes weaker when 
the flow roughness is increased. This prediction is con- 
firmed by the numerical studies presented in next Sec- 
tion. 

Finally it should be pointed out that although the 
change of variables © - ([7]) can be applied equally in 
three dimensions, the above analysis does not carry over 
to higher dimensions. Firstly, as already pointed out, 
an additional drift term arises. This Ito-term renders 
a straightforward derivation of an analytic solution for 
the deterministic drift impossible. Secondly, for higher 
dimensions the fixed point of the reduced dynamics is 
located far from the origin, see p^ . Hence the approxi- 
mations made above for d = 2 are not applicable. Careful 
numerical studies are needed to understand whether or 
not algebraic tails are also present in higher dimensions. 

III. CORRELATION DIMENSION AND 
APPROACHING RATE 

Particle clustering is often quantified by the radial dis- 
tribution function g{r), which is defined as the ratio be- 
tween the number of particles inside a thin shell of ra- 
dius r centered on a given particle and the number which 
would be in this shell if the particles were uniformly dis- 
tributed. This quantity enters models for the collision 
kernel [121]. Following [S, [H, [H, [H] , we consider a dif- 
ferent, but related way to characterize particle clustering. 
Instead of the radial distribution function we evaluate 
the correlation dimension I?2 of the set formed by the 
particles. This dimension is widely used in dissipative 
dynamical system theory and in fractal geometry (see, 
e.g., [H, [i^])- It is defined as the exponent of the power- 
law behavior at small scales of the probability P2{r) of 
finding two particles at a distance R<r: 

I?2 = limdsW, d2(r) = ^^J^^, (13) 
r^o d in r 

where the logarithmic derivative d2{r) is called the local 
correlation dimension. I?2 relates to the radial distribu- 
tion function via Ing(r)/ Inr — > 7)2— d for r—^0. For uni- 
formly distributed particles, T>2—d, so that g(r) = 0(l). 
On the contrary, when particles cluster on a fractal set, 
7)2 <d and q(r\ diverges for r^O. This was also found 
numerically in [17| . 

Depending on whether the carrier flow is spatially 
smooth {h = 1) or rough {h < 1), I?2 and d2{r) behave 



differently. In the former case, random dynamical sys- 
tem theory [2l[ suggests that within the 2 x d position- 
velocity phase space, particles converge onto a multifrac- 
tal set with correlation dimension < I?2 < 2d. Here 
I?2 denotes the correlation dimension in the full phase 
space. It is defined in complete analogy to I?2 through 
the scaling behavior of the probability P2{'r) to find two 
particles at a distance less than r in phase space: 

P2(r)~r'^= for r^O. (14) 

The distance r is now computed by using the phase-space 
Euclidean norm + |V/Z3ip; V is normalized by 

the typical fiuid velocity gradient Di for dimensional rea- 
sons. The physical-space correlation dimension I?2 is ac- 
tually the dimension of the projection of the set from 
the full phase space onto the position space, and it is 
also expected to be fractal (see Section I VIII for details 
on the relation between I?2 and 7)2) ■ We focus in this 
Section on quantifying clustering in position space and 
hence consider only 7)2 and d2{r). 

Balkovsky et al. argued in [43| that particles do not 
form fractal sets in non-smooth flows because the cor- 
relation function of the particle density fleld should be 
a stretched exponential. Clustering and inhomogcneities 
are hence not quantified by a fractal dimension but by the 
detailed scale dependence of d2{r). However, as discussed 
in the Introduction, one expects the statistical properties 
of two particles separated by a distance r in a flow with 
Holder exponent h to depend on the local Stokes number 
St{r) = DiT /r'^'^^~^'^ only, which for smooth flows degen- 
erates to a scale independent number, St(r) = St = Dit. 
In rough flows, at scales small enough, particles move bal- 
listically and distribute homogeneously as the Lagrangian 
motion is too fast for the particles to follow (St(r)^oo 
as r ^ 0) and hence I?2 = d for all particle response times 
T. However, information on the inhomogcneities of the 
particle distribution at larger scales can still be obtained 
through the scale-dependence of the local correlation di- 
mension d2{r) defined in (|13p . 

The relevance of the local Stokes number and of the 
local correlation dimension is confirmed by numerical ex- 
periments of planar suspensions. Simulations were per- 
formed by directly integrating the reduced system de- 
scribed in previous Section. Figure [3] shows d2 {r) as a 
function of St(r) for various values of h. The curves ob- 
tained with different values of the response time r col- 
lapse onto the same /i-dependent master curve once the 
scale dependency is reabsorbed by using St(r). In the 
plot, only scales far from the boundaries were considered, 
as otherwise the self-similarity of the fiuid fiow is broken. 
The data ioY h — \ estimate the limit of d2{r) as r— >0, 
and so correspond to the value of the correlation dimen- 
sion 7)2- As anticipated in the previous Section, Fig. [3] 
also shows that clustering is weakening when the rough- 
ness of the fluid velocity increases (i.e. when h decreases). 
In particular, minr{(i2(^)} gets closer to d, i.e. particles 
approach the uniform distribution as /i ^ 0. Finally no- 
tice that for St(r) — > 0, i.e. at large scales in rough flows, 
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d2 (r) d as well. This is due to the fact that at these 
scales the Lagrangian motion becomes much slower than 
the relaxation time of the particles. The particles thus re- 
cover the tracer limit and distribute homogeneously. As 
we will see in Section IVll the local dimension d2{r) tends 
linearly to the space dimension d when St(r) —>■ with 
a factor whose dependence on h and d can be obtained 
analytically by perturbative methods. 

The radial distribution function and hence the correla- 
tion dimension give only partial information on the rate 
at which particles collide. Indeed, in order to evaluate 
the collision rate, one needs to know not only the proba- 
bility that the particles are close to each other, but also 
their typical velocity difference. Here, following we 
study the approaching rate ^(r) defined as the flux of 
particles that are separated by a distance less than r and 
approach each other, i.e. 

K{r) = {R ■ R/\R\e{-R ■ R/\R\) e{r - \R\)) , (15) 

where 8 denotes the Heaviside function and the average 
is defined on the Lagrangian trajectories. As detailed 
in ^1^, K{r) is related to the binary collision rate in the 
framework of the so-called ghost collision scheme [23| . 
Within this approach collision events are counted while 
allowing particles to overlap instead of scattering. At 
small separations, ^(r) behaves as a power law. This 
algebraic behavior allows defining a local Holder exponent 
7(r) for the particle velocities 

7W = ^-rf2(r). (16) 
inr 

In the above definition the contribution from clustering, 
accounted for by the local correlation dimension d2{r), 
is removed. The local Holder exponent 7(r), similarly 
to d2{r), tends to a finite hmit F as r ^ which, for 
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FIG. 3: Local correlation dimension d2{r) versus the scale- 
dependent Stokes number St(r) — Dir/r^'^"''' for two- 
dimensional flows with different h. Symbols denote different 
particle response times r. For h — 1, 1^2 = d2(r 0) is 
displayed and St(r) = St = Dir. 



particles suspended in a smooth flow {h = 1), depends 
non-trivially on the Stokes number. 




St(r) 

FIG. 4: Ratio between the local Holder exponent 7(r) of the 
particle velocity and that of the fluid h versus St(r). The 
symbols in each curve refer to different values of the particle 
response time r. As in Fig. [S] for h = 1, the small scale 
limiting value F is depicted. 

Figure [3] shows numerical estimations of j{r)/h as a 
function of St(r) for various values of h. In the smooth 
case {h = 1), the limit value F decreases from F = 1 for 
St = 0, which corresponds to a differentiable particle ve- 
locity field, to F = for St^oo, which means that parti- 
cles move with uncorrelated velocities [16| . The fact that 
F < 1 is due to the contribution of caustics appearing in 
the particle velocity field 0, El [13, [H, HI (see Sect.|V] 
for a discussion in c? = 1). Similarly, in non-smooth flows 
7(r) is asymptotically equal to the fluid Holder expo- 
nent h at large scales (St(r) — > 0), and approaches at 
very small scales (St(r)^oo). Therefore, all the relevant 
information is entailed in the intermediate behavior of 
7(r). The latter should only depend on the fluid Holder 
exponent and on the local Stokes number, as confirmed 
by the collapse observed in Fig. HI Note that the tran- 
sition from 7(r) = h to 7(r) = shifts towards larger 
values of the local Stokes number and broadens as h de- 
creases. The fact that 7(r) = h ioi r ^ oo implies that 
the particles should asymptotically experience Richard- 
son diffusion just as tracers (see Sect.|lV]for details). For 
comments on how the findings reported in this Section 
translate to realistic turbulent flows, we refer the reader 
to Section lVlIIl 



IV. STRETCHING RATE AND RELATIVE 
DISPERSION 

This Section is devoted to the study of the behavior of 
the distance R{t) between two particles at intermediate 
times t such that i?(0) <C R{t) <C L. For convenience, 
we drop the reflective boundary condition at R — L and 
consider particles evolving in an unbounded domain. 
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FIG. 5: Lyapunov exponent versus St: the circles are the 
numerical measurements while the dashed line corresponds 
to Eq. (|20|) . Inset: rate function H associated to the large 
deviations of the stretching rate for three values of St; the 
solid line corresponds to H for tracers for, whose analytic 
expression is known (see, e.g., [3]). 



We first consider a difTerentiable fluid velocity field 
{h = 1). In this case, the time evolution of the distance 
R{t) is given by so that 



R{t) = i?(0) exp 



cri(i')dt' 



(17) 



and the particle separation can be measured by the 
stretching rate ^{t) = {l/t)\n[R{t)/R{0)]. It is assumed 
that the reduced dynamics ^ - PH)) is ergodic. There is 
currently no rigorous proof of ergodicity. However, such 
an assumption relies on numerical evidence and on the 
following phenomenological argument. The deterministic 
loops described in Section |TT] are randomly initiated by 
the near-origin behavior of the system, providing a mech- 
anism of rapid memory loss that might ensure ergodicity. 
With this assumption, the time averages converge to en- 
semble averages, so that 



ai{t')dt' 



^1 



as t 



(18) 



In other words, the distance between particles asymptot- 
ically behaves as R{t) = R{0) exp{tX), where A = (cti) is 
a non-random quantity referred to as the Lyapunov ex- 
ponent. A positive Lyapunov exponent implies that the 
particle dynamics is chaotic [l9| . 

Figure O shows numerical measurements of the Lya- 
punov exponent A. The exponent remains positive for all 
values of the Stokes number. This means in particular 
that particles suspended in incompressible flow cannot 
experience strong clustering, which consists in the con- 
vergence of all trajectories together to form point clus- 
ters. This contrasts with the case of compressible flows 
where, for suitable values of St and of the compressibil- 
ity, negative Lyapunov exponents are observed A 
first attempt to derive an analytic expression for A (St) 
was proposed by Piterbarg [l3|. His approach is based 



on studying the Laplace transform Lp{p) of the distribu- 
tion of the complex random variable z — ai + ia2, i.e. 
tp{p,t) = (exp(— pz(i))) which satisfies 



dtip^ ~{p/t) dp^+pd^pip - (2i?i/T)p2<^. 



(19) 



If ip{p, t) reaches a steady state at large times, one can 
infer an analytic expression for the asymptotic solution 
Voo(p) by requiring that the right-hand side of (|19p van- 
ishes. It is then straightforward to deduce that the Lya- 



punov exponent satisfies A 
implies 



-limp^o 5?{c?p</'oo}- This 



A 



El 

2St 



m 1 



Ai'(x) 



(16 St) 



-2/3 



(20) 



^/x A{x) 

where Ai and Ai' designate the Airy function of the first 
kind and its derivative respectively. This prediction is 
compared to the numerical measurements in Fig. [51 As 
stressed in [IBj, there is evidence that the moments Lp{p, t) 
do not converge to a steady state, but rather diverge at 
large times. This might explain the discrepancies ob- 
served in Fig. [5l However, the numerical precision is not 
high enough to test the presence of corrections to the 
analytic expression 

At large but finite time i, the distance between the 
two particles is measured by the stretching rate ^(t) = 
{l/t)ln[R{t)/R{0)]. This quantity becomes more and 
more sharply distributed around the Lyapunov exponent 
A as t increases. More precisely, it obeys a large devi- 
ation principle and its pdf p(/i, t) takes the asymptotic 
form (see, e.g., [1|) 

jliip{ti,t)^-H{p), (21) 

where _ff is a positive convex function attaining its min- 
imum in /i = A, in particular H(X) = 0. The rate func- 
tion H measures the large fluctuations of /i, which are 
important to quantify particle clustering. Rate functions 
obtained from numerical experiments are represented in 
Fig.[5]for various values of the Stokes number. The func- 
tion becomes less and less broad when St increases, a phe- 
nomenon that can be quantified in the limit St — s- oo as 
discussed in Section IVlII Note that the same qualitative 
behavior is also observed for heavy particles suspended 
in homogeneous isotropic flow 
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We now turn to the case of particles suspended in non- 
differentiable flows {h < 1). As we dropped the bound- 
ary condition, the initial inter-particle distance R{0) is 
the only relevant length scale. By using i?(0) instead of 
L in the change of variables (O - ([7]) the problem of rel- 
ative dispersion is expressed solely in terms of the the 
Holder exponent h and of a time-dependent Stokes num- 
ber which can be defined in terms of the local Stokes 
number as Stt = DiT/[R{t)]'^^^^'^\ In particular, the 
evolution of R{t) directly follows from the initial its value 
Stg. From the evolution equation (fTU)) for the reduced 
separation p{t) = [R{t) / R{0)]^~'^ , we obtain 



p{t)^l + {l-h) / ai{t')dt', 



(22) 
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FIG. 6: Time evolution of the average rescaled separation 
{{p{t) ~ p(0))) for different initial Stokes numbers Sto, and 
h = 0.4, 0.6, and 0.8 (from top to bottom). Inset: long- 
time behavior of the time-dependent Stokes number Stt — 
D\T/p^{t) for different Sto and the same three values of h 
(now from bottom to the top). The segments on the left 
indicate the slopes —2/(1 — h) corresponding to the regime of 
ballistic separation. 



where p(G [0, oo)) typically increases with time. The time- 
dependent Stokes number St* = DiT/i?2(i-'») ^ Sto/p^, 
which measures the effect of inertia when the particles 
are at a distance i?(i), decreases with time. Hence, con- 
versely to the case of differentiable carrier flow, ai is not 
a stationary process and the integral in (|^^ does not 
tend to t{(7i). 

Hereafter, we confine the discussion to the case Sto >• 1 
because it contains a richer physics than smaller Sto. As 
observed from Fig. [SI we can distinguish two regimes in 
the time behavior of p{t). At first the particle separation 
evolves ballistically, i.e. R{t) (x t, meaning that the time- 
dependent Stokes number Stt decreases as (see 
inset of Fig. and reaches order-unity values for t k t. 
During this phase, the time growth of p is accelerated 
or slowed down and ultimately reaches a diffusive be- 
havior (X t^^'^. This corresponds to the limit of tracers, 
which is approached when Stt «C 1. At this stage, the 
inter-particle distance behaves as R{t) oc and, 
consequently, the Stt decreases as 1/t (see Fig. |6l). 

The convergence to tracer diffusion in the limit of large 
distances R gives an original way to interpret Richard- 
son's law for delta-correlated velocity fields in terms of 
the asymptotic behavior of the reduced variables (H)) - (O . 
When p is large, the quadratic terms in the drift of equa- 
tion ([8]) can be neglected and ci behaves as an Ornstein- 
Uhlenbeck process with response time r. However, when 
(Ti becomes of the order of p/(/it), the quadratic terms 
cease to be negligible and they push the trajectory back 
to ai > 0. This process happens on time scales that 
are of the order of unity and thus much smaller than the 
time scales relevant for large-scale dispersion. Hence the 
dynamics of ai{t) can be approximated as an Ornstein- 



FIG. 7; Pdf of the rescaled separation p{t) for various combi- 
nations of Sto and large times t. The solid lines represent the 
limiting distribution given by (|23p with A = 1/4. 

Uhlenbeck process with reflective boundary condition on 
(Ti — p/{hT). This implies that p has a diffusive behav- 
ior. More speciflcally, numerical simulations (see Fig. [7]) 
show that the pdf of p behaves as 

p(p,i)cxp''<-(''+i)/2exp[-VA], (23) 

where v = {l + h)/{l — h) and yl is a positive constant. At 
large times and consequently large distances Stt 0, the 
tracer limit is fully recovered as conflrmed by expressing 
the above relation in terms of the physical distance R = 
pi/{i-h)^ jn(jgg(j it becomes identical to the law that 
governs the separation of tracers in a Kraichnan flow [28j . 
However, a direct derivation of ([25]) in terms of the p and 
(T dynamics is still lacking. 

V. EXACT RESULTS IN ONE DIMENSION 

A number of analytical results were derived for one- 
dimensional flows [1^, [13, dll . Although such flows are 
always compressible, their study helps improving the in- 
tuition for the dynamics of inertial particles in higher- 
dimensional random flows. In particular, several results 
on caustic formation hold also in two-dimensional (in- 
compressible) flows because the typical velocity fluctu- 
ations, which lead to caustic formation, are effectively 
one-dimensional. 

Here, we focus on one-dimensional smooth flows, for 
which the equations analogous to ([5]) - PU)) reduce to 

& = -cr/T-cT^ + VCri(t), (24) 
R = aR, (25) 

where a = V/R and, as in (Ell-dTOl), C = 2Di/t^. 
The quadratic term in ([M]) implies that a can escape to 
— oo with a finite probability. These events are the one- 
dimensional counterpart of the loops described in Sec- 
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tionlTTland correspond to the formation of caustics: par- 
ticle trajectories intersect with a finite relative velocity. 
Note that the equation for a decouples from the equation 
for R, so that it can be studied separately. Stationary 
statistics of a can be described by the pdf P{a) which 
obeys the one-dimensional Fokker-Planck equation 

[d,{cTlT + a^)+{Cl2)dl]P{cj) = Q. (26) 

This equation can be rewritten as da-J{cr) — 0, where 
J(cr) = (cr/r + cr2)P((7) + CP'{a)/2 is a probability flux 
in the u-space. Equation ([26]) is supplied by the bound- 
ary conditions J(-l-oo) = J{—oo), which are required 
to resolve escapes to infinity and thus caustic forma- 
tions. Indeed such events correspond to particle cross- 
ings during which R and V remains finite, so that 
a — V/R changes sign. Hence, all particles escaping to 
(7 = -|-oo reappear at a — —oo. The stationary solu- 
tions of Eq. ([^ satisfying such a boundary condition 
corresponds to a constant flux J and can be written as 

P(^a)^Me-^uia)/c fdaV^^^')/^, (27) 

^ J-oo 

where U{a) — -I- (t^/2t. Note that as in two dimen- 
sions, P{<j) has power-law tails. The argument presented 
in Section [III can actually be straightforwardly applied 
with the difference that there is no loop anymore but just 
escapes to infinity occurring with a probability that is in- 
dependent of a. This leads to P{a) cx \<t\^^ for |cr| — > oo 
(the exponent is actually —{l + l/h) in the general case 
of Holder-continuous carrier flows). 

Using the constant-flux solution P?)) . one can de- 
rive the Lyapunov exponent A — (a). As shown in 
(30| . its value non-trivially depends on the Stokes num- 
ber. For St = DiT <C 1, it is negative and behaves 
like A ~ — -Di while for St 1 it becomes posi- 
tive and its value is given by the asymptotic expression 
A ~ i:)iSt-2/3V3125/6r(5/6)/(24V^) > . There ex- 
ists a critical value of the Stokes number (« 0.827) for 
which the Lyapunov exponent changes its sign. This 
phenomenon of sign-changing has been dubbed path co- 
alescence transition by Wilkinson and Mehlig in [30| . It 
is closely related to the aggregation-disorder transition 
discussed in [2^. The sign of the Lyapunov exponent 
determines how the distance between two initially close 
particles evolves with time. It turns out that the answer 
depends on the particle size: small particles (with small- 
enough Stokes numbers) tend to approach each other, 
while large particles (with large Stokes numbers) get sep- 
arated by the flow. 

Another important phenomenon which was extensively 
studied within the one-dimcnsional model is the for- 
mation of caustics. The average rate of caustics for- 
mation is given by the absolute value of the probabil- 
ity flux J. For large values of the Stokes number it 
can be written as |J| ~ DiSt-'^^^T{5/6)12^^^/{8Tr^^^), 
while for small Stokes it becomes exponentially small 



I J| - Di(27rSt)-i exp[-l/(6St)]. The formation of caus- 
tics is a stochastic process, whose properties can be de- 
scribed by the pdf of the caustic formation time T. In 
(Slj it is shown that for St <C 1 this pdf can be estimated 
as P(r) cx exp[-l/(6St)] for r < T < r exp[l/(6St)] 
and P{T) cx exp [~w/{3CT^)], with w = r(l/4)8/967r2 
(r denoting the Gamma- function here), for T <^t. The 
exponential factor exp[— l/(6St)] which characterizes the 
small rate of caustic formations for St <C 1 can be easily 
explained if one formally considers Eq. ([M)) as a Langevin 
equation for a particle which is driven by the thermal 
noise •q{t) and evolves in the potential U{<t). In this 
case, the rate of caustic formation is given by the prob- 
ability for the particle to tunnel through the potential 
barrier in U (cr) . Such probability can be estimated as 
exp[— l/(6St)]. For large Stokes numbers, the barrier 
disappears and the rate of caustic formation is not ex- 
ponentially damped anymore. 



VI. SMALL STOKES NUMBER ASYMPTOTICS 

This Section reports some asymptotic results related 
to the limit of small particle inertia. The flrst part sum- 
marizes the approach developed by Mehlig, Wilkinson, 
and collaborators for differentiable flows (/i= 1). In anal- 
ogy to the WKB approximation in quantum mechanics 
(see, e.g. [s^]), the authors construct perturbatively the 
steady solution to the Fokker-Planck equation associated 
to the reduced system ([S])-®. In the second part of 
this Section original results are reported where the par- 
ticle dynamics is approximated as the advection by a 
synthetic flow comprising an effective compressible drift 
which accounts for leading-order corrections due to par- 
ticle inertia. 

Mehlig and Wilkinson proposed in [23| (see also [33| ) 
to approach the limit of small Stokes numbers in terms 
of the variables xi — (j j D-i)^l'^a\ and = {t /'iDiY^'^(j2- 
From equations®-©, their time evolution follows to 
satisfy 

= -xi-£[xl~2,xf^+V2-qi{s), (28) 
±2 = -X2-2exiX2 + V2ri2{s) , (29) 

where e = -s/St, dots denote derivatives with respect to 
the rescaled time s = t/T and rji and 772 are independent 
white noises. The evolution equations ((28| - ((29)) can be 
written in vectorial form, namely x = —x-\-eV{x)+^/2r|^ 
where x = {xi,X2), t) = (771,772) and V denotes the 
quadratic drift. The steady-state probability density 
p{x) is a solution to the stationary Fokker-Planck equa- 
tion 

VIp + V, • {xp) = eV. • [\{x)p\ . (30) 

Next step consists in writing perturbatively the proba- 
bihty density of x as p[x) = exp(— |a;p/4) (Qo + eQi + 
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^2 • • • )• The functions Qfc satisfy the recursion rela- 
tion HoQk+i = HiQk, where 



Ho 



l + V^ -|a;|V4, 
V^-V(a;)+a;-V(a;)/2. 



(31) 
(32) 



The operator TLq is the Hamiltonian of an isotropic two- 
dimensional quantum harmonic oscillator. This suggests 
introducing creation and annihilation operators and to 
expand the functions Qk in terms of the eigenstates of 
the harmonic oscillator (see [13, [1^ for details) . 

This approach yields a perturbative expansion of the 
Lyapunov exponent [2^ 



A 



Di{xi)/e = 2Di J2 a^e'^' = J2 "fe^*^' (^3) 



fc>0 



fc>0 



where the coefficients ak satisfy the recurrence relation 



flfc+i = 4(4 - 3k)ak - 2 ^ 



(34) 



£=0 



with Co = 1. For large fc, these coefficients behave as 
flfe ^ (— 12)'^fc!, so that the series ((33)) diverges no matter 
how small the value of e (and thus of St). Hence the sum 
representation of A makes sense as an approximation only 
if truncated at an index fc* for which |afcSt'^| attains its 
minimum. For small values of St, ki, l/(12St) and the 
error of the asymptotic approximation is of the order of 
the smallest term, namely ~ jafe^St'"'* | ^ exp[— l/(12St)]. 
This approach was refined by Wilkinson et al. [33| adopt- 
ing an approach based on Pade-Borel summation, which 
was found to yield satisfactory results. 

The non-analyticity of A(St) at St = is interpreted in 
[131 as a drawback of the perturbative approach. Indeed 
the quadratic terms in ([28|l - ([29]) are not negligible for 
all values of xi and X2'- When \x\ becomes larger than 
e^^ they are actually dominant and the trajectory per- 
forms a loop in the x (or cr) plane (see Section|lIl). When 
St = is small, the probability to initiate such a loop 
is given by the tail of the distribution governing scales 
\x\ <C e~^, and is hence cx exp[— 1/(6£^)], which coin- 
cides with the one-dimensional result discussed in pre- 
vious Section, confirming the relevance of d = 1 physics 
to the formation of caustics in higher dimension. Tak- 
ing into account this correction due to caustics, i.e. the 
contribution of events when the particles approach very 
close to each other keeping a finite velocity difference, 
Mehlig and Wilkinson proposed to write the Lyapunov 
exponent as 



-l/(6St) 



2 UkSt'^ 

k=0 



(35) 



where S is a positive constant. We finish this summary 
by stressing that this approach equally applies to the case 
of compressible carrier flows [2J|, and was extended to 
three dimensions where it yields a prediction on the St- 
dependence of the three largest Lyapunov exponents [33j . 



The above perturbative approach can be generalized to 
small particles evolving in rough flows. For small (local) 
Stokes numbers, the characteristic time scales of veloc- 
ity evolution are much smaller compared to the temporal 
scales associated to the dynamics of the particle separa- 
tion. Therefore, one can obtain the effective equation for 
the evolution of particle separation by averaging over the 
fast velocity difference variables. The systematic mathe- 
matical strategy of such an averaging was proposed in [s^l 
in the context of stochastic climate models. This strat- 
egy is closely related to the Nakajima-Zwanzig technique 
which was developed to study similar problems arising in 
damping theory |35l . [36j . Applications of this technique 
to the elimination of fast variables in Fokker-Planck equa- 
tions are discussed in , '38^1 . In this framework one can 
derive an expansion for the Fokker-Planck type operator 
entering into the equation for the slow-variable proba- 
bility distribution function. In our case, this leads to a 
closed equation for the pdf of the particle separation R. 
This equation can be used to determine the local corre- 
lation dimension d2{r) for St(r) <C 1. We present here 
only the general idea and the main results; details of the 
calculations will be reported elsewhere. 

To carry out the above-mentioned procedure the joint 
position- velocity pdf p{r, v) is approximated by 



p{r, v) ~ p{r)Pr{v) + p{r, v), 



(36) 



where p{r, v) denotes subleading terms which are 0(St); 
Pr{v) is the stationary distribution associated to the fast 
velocity variables and satisfies the Fokker-Planck equa- 
tion 



LoPriv) 



1 ■ . b'-Hr) ■ ■ 

' 9 V 



Priv) = 0, (37) 



with the normalization condition jdv Pr{v) = 1. With- 
out loss of generality, it is assumed that the subleading 
terms p(r, v) in the approximation (j36p do not contribute 
to the normalization condition, so that Jdvp{r,v) ~ 
p{r). The effective equation for p{r) can be derived by 

introducing the expansion p(r) = J2T=o St'^''^Pfe('")- This 
expansion, which enters the definition p6p . is then sub- 
stituted into ([4]) and all terms of the same order in St 
are collected. Note, that the operator Li = d^v^ entering 
Eq. (HI) is smaller than the other operators by a factor 
St^/^. The chain of equations for Pfe(r) has a solvability 
condition that results in the following effective equation 
for p{r): 



(38) 



Afi + Af2 + • • • pir) = 0, 



where the operators Mk can be written as 



MkPir)^ dv (LiLo^) Lip{r)Pr{v). (39) 



Lq ^ denotes here the inverse of Lq, i.e. the Green func- 
tion obtained from (|37p with the right-hand side replaced 
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by a 5 function. This operator is defined in such a way 
that JdvL^^f{v) = for any function f{v) satisfying 
Jdv f[v) = 0. One can check that the leading-order oper- 
ator is Ml = dlV^{r)dl which, as expected, corresponds 
to turbulent diffusion. Indeed the dynamics of tracers is 
recovered when St 0. The pdf p{r) which solves the 
equation Mip{r) = is simply the uniform distribution. 
To measure particle clustering, which can be estimated 
for instance by the local correlation dimension d2 (r) (see 
Section IIIII) , one has to calculate the next order opera- 
tors. It can be easily checked that all operators Mk of 
even order k are zero. The first non- vanishing correc- 
tion to Ml is thus given by the third order operator Afa. 
When interested in the stationary distribution only, the 
terms which enter this operator and which are associated 
to transients can be disregarded and one can write: 

Afg . = dl[V'-l with V' = -'^ {d^dlb'^) {d^b^^) . (40) 

The operator M^ can be interpreted as an effective drift 
in r-space and, for the Kraichnan model, represented as 
V' ^ -2{<f - l){d-2 + 4:h)h^St^{ry. The functional 
form of this drift implies that the first non- vanishing cor- 
rections to the uniform distribution are proportional to 
St(r). Indeed, for isotropic flows one can look for a so- 
lution p{r), which depends only on the modulus r of 
its argument. In this case Eq. ([55]) becomes an ordi- 
nary differential equation of Fokkcr-Planck type. Look- 
ing for a non-flux solution one readily obtains the de- 
sired p{r). In rough flows {h < 1), one has Inp(r) ~ 
[{d +l){d-2 + 2h)h'^/{l - h)] St(r) and the local corre- 
lation dimension behaves as 

.,(.).. -Ml±l)(l^l±i^St(.). (41) 
^ ' d-2 + 2h ^ ' ^ ' 

Note that the second term on the right-hand side of the 
above expression disappears for h 0, confirming once 
again the finding of the previous Sections about the de- 
crease of clustering going from smooth to rough flows. 
For differentiable carrier flows {h ~ 1), the distribution 
has algebraic tails: lnp{r)^ —2{d + l){d + 2)Stlnr, and 
hence the correlation dimension behaves as 

I?2 = d-2(d+l)(rf + 2)St + 0(St2). (42) 

The dimension deficit d — I?2 is equal to 24St for two- 
dimensional flows and to d — I?2 — 40St for three- 
dimensional ones. The latter result is in agreement with 
the dimension deficit of the Lyapunov dimension reported 
by Wilkinson et al. in [33| . The above predictions on 
the dimension deficit, for smooth fiows, are in very good 
agreement with numerical simulations in two and three 
dimensions, see Fig. HI We conclude this Section by 
noticing that in time-correlated random smooth flows, 
as well as in developed turbulence, the dimension deficit 
has been shown to be oc St^ 0, [III, [H, S^l . Therefore, 
including temporal correlations seems to be crucial to re- 
produce the details of the small-Stokes statistics of tur- 
bulent suspensions. 
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FIG. 8: Dimensional deficit 2 — D2 versus St in d = 2 for 
smooth flows {h = 1). Inset: same for d = 3. Points rep- 
resent numerical results and the straight line corresponds to 
the perturbative predictions given by (|42|l for d = 2 and 3 
respectively. 



VII. LARGE STOKES NUMBER 
ASYMPTOTICS 

Particles with huge inertia (St ^ 1) take an infinite 
time to relax to the velocity of the carrier fiuid. They 
become therefore uncorrelated with the underlying flow 
and evolve with ballistic dynamics, moving freely and 
maintaining, almost unchanged, their initial velocities. 
This limit is particularly appealing for deriving asymp- 
totic theories [l^. In this Section, we focus on two as- 
pects, namely the problem of the recovery of homoge- 
neous/uniform distribution for St ^ 1 and the problem 
of the asymptotic scaling for the statistics of the particle 
separation and of the velocity differences. 

A. Saturation of the correlation dimension 

Ballistic particles injected homogeneously and uni- 
formly remain so [4l[. Hence for the correlation di- 
mension associated with their distribution \Vi\ one has 
I?2 = d. This result follows directly from the Fokker- 
Planck equation (|4]), which can be seen as an advection- 
diffusion equation in phase space. The effective flow is 
compressible because of the term —d^vjT but, in the 
limit S't — !■ cx), it becomes negligible and the equation 
reduces to diffusion plus advection by an incompressible 
flow. The resulting stationary pdf is thus uniform in 
phase space and hence in its projection in position space. 
Moreover, as particle velocities and fluid flow are uncor- 
related and consequently the particles are not correlated 
with each other, the exponent F which characterizes the 
small-scale behavior of the approaching rate (see Sec- 
tion |TTT| vanishes. Thus D-i^d and F^O for St ^00. 

This asymptotic regime can be achieved via two possi- 
ble scenarios: (a) asymptotic convergence of I?2 to d, and 
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(b) saturation of ©2 to d for Stokes numbers above a crit- 
ical value St^. In what follows, we provide evidence for 
(b), limiting the discussion to two-dimensional smooth 
flows. 

Let us first discuss a phenomenological argument in 
favor of saturation. As already noted in Section HIIl their 
dissipative dynamics yields the phase-space trajectories 
of the particles to converge onto a random, dynamically 
evolving attractor, which is typically characterized by a 
multifractal measure [l^, [lO]. In our setting, this mea- 
sure is the phase-space correlation dimension defined in 
equation . Ballistic motion for St ^ 1 corresponds 
to 2?2 2d, therefore a critical Stokes number St^ exists 
such that 2?2(St^) = d. The particles' spatial distribution 
is obtained by projecting the (2 x c?)-dimensional phase 
space onto the d-dimensional physical space. It is tempt- 
ing to apply a rigorous result on the projection of random 
fractal sets [22','42'| stating that for almost all projections, 
the correlation dimension of the projected set is related 
to that of the unprojected one via the relation 



I?2 = min{(i, 7^2} ■ 



(43) 



Having I?2(St^) = d with the above expression implies 
that I?2(St) = d for all St > St^. Unfortunately, there is 
a priori no reason for assuming some kind of isotropy in 
phase space which justifies the validity of (|43| . We thus 
proceed numerically. 

As Eq. (|43|) requires the isotropy of the set, we have 
tested whether this applies to our case. The correla- 
tion dimension of different two-dimensional projections 
was evaluated through the computation of the probabil- 
ities -P^''^(r) of having two particles at a distance less 
than r using the norm ^ — (5^ + (5^, with a,/3 = 
X,Y,Vx / Di,Vy / Di, and 5a denoting the coordinate-a 
separation between the two particles. Note that a = X 
and P — Y corresponds to the spatial correlation di- 
mension discussed so far. Figure [5] shows the logarith- 
mic derivatives (dlnP2"'^(''))/('ilnr) for various a, /3 and 
three different values St. All curves collapse within error- 
bars, confirming that the projection is rather typical and 
thus strengthening the argument in favor of saturation. 
However, as can be seen in Fig. [HI the logarithmic deriva- 
tives on the different projections are curved, indicating 
behaviors different from the expected power law. It is 
therefore difficult to decide whether or not the satura- 
tion occurs. As discussed in [3], one can understand 
the curvature of the local slopes with the presence of 
sub-dominant terms, e.g., with the superposition of two 
power laws P2{r) ~ Ar°' + Br^. In our case, one can 
expect that 



P2(r) = Ar^^ +Br'^ 



(44) 



where d and 2?2 are the only dimensions entering the 
problem [l^. For I?2 < d, the second power law can 
be interpret also as the contribution of caustics [H, |2^ : 
With non-zero probability, particles may be very close to 
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FIG. 9: Logarithmic derivative (dlnP^ (r))/(dlnr) for dif- 
ferent projections a, 13 for St = 0.5, St = 1 (shifted up by a 
factor 1), and St = 1.5 (shifted up by a factor 2). A small mis- 
match in the scaling range can observed for large r (this is un- 
avoidable as positions and velocities involve different scales). 



each other with quite different velocities, see SectionlVl 
Once projected onto physical space, caustics appear as 
spots of uncorrelated particles, and hence, the correlation 
dimension is locally I?2 = d. The validity of (|44p as 
well as of the projection formula ([43]) was confirmed in 
Ref. [H, . 
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FIG. 10: Physical space ©2, and phase-space ©2 correlation 
dimensions versus St as obtained by using (|44|l for fitting the 
exponents. Errors are of the order of the size of the symbol. 
The arrow indicates the estimated location of St^ . 

Figure [TO] summarizes the results depicted above. In 
particular, I?2 clearly displays a crossover to values larger 
than d for St > St^ w 0.6. I?2, once properly estimated 
by using (|43p . displays the saturation to d = 2 above St^, 
at which the large Stokes asymptotics starts, at least for 
the particle distribution. 

Let us comment briefly on the implication of satura- 
tion on the behavior of the approaching rate which, in the 
limit St —> 00, is characterized by the exponent F ^ 0. 
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Similarly to 'D2^ deviations of V from its limiting value 
cannot be determined by scaling arguments. Saturation 
of would however affect F. This is related to the dom- 
inant contribution of caustics which might imply also the 
saturation of F to for sufficiently large Stokes numbers. 
Though numerical experiments confirm this scenario , 
saturation cannot be studied with as much detail as for 
I?2- At present, there is no simple phenomenological ar- 
gument for the subleading terms as for I?2- 



B. Scaling arguments 

The limit of large values of the Stokes number can 
be approached by assuming r — > cxd and keeping C — 
2Di/(tL1~^)2 constant. The dynamics dH) - Q for the 
relative velocity differences can then be approximated by 

di ~ - {hal - al) /p+^/Cm, (45) 
&2 ~ -(/i+l)aiCT2/p+ V(l + 2/i)C772 . (46) 

For a given exponent /i, the limiting dynamics depends 
solely on C while — after non-dimensionalizing time and 
relative velocities by r — the general dynamics depends 
on St(L) only (see the Introduction). This congruence, 
which was first used in [45J for determining the large-St 
behavior of the Lyapunov exponent, allows to derive scal- 
ing arguments of various other quantities characterizing 
two-particle dynamics. 

Let us detail this for the distribution of the longitu- 
dinal velocity difference ai. It is clear from the above 
considerations that for fixed h and cti ^ (I/t) the fol- 
lowing relation holds 

rp(Tai;St) ~p(ai;C). (47) 

Differentiating with respect to Di and r gives a necessary 
condition for such a behavior: p must satisfy 

p + Gid^.p + iCdcP^G, (48) 

which itself implies p(cti; C) = C-^'^ f{C-^'^ai), so that 

p{ai) ~ St-^/V/(St-^/Vcri) for St > 1. (49) 

As shown in Fig. [Tl] this asymptotic scaling behavior 
can be observed numerically. As a consequence of (|49|) . 
for differentiable carrier flows {h — \) the Lyapunov ex- 
ponent A = ((Tl), which measures the asymptotic growth 
rate of the inter-particle distance (see Section lIVp , be- 
haves as 
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FIG. 11: Pdf of the non-dimensional longitudinal velocity dif- 
ference (Tl at large values St (symbols are for different values) 
for various values of h. 

is shown in [l^ that this result also holds in three di- 
mensions. Its confirmation by numerical simulations is 
illustrated in Fig. [12] 

The scaling argument described above can be car- 
ried forward to the fluctuations of the stretching rate 
^{t) = (l/i)ln[i?(t)/i?(0)]. As we have seen in Sec- 
tion [iVl for large times the distribution of ^ obeys the 
large deviation principle (|2ip . It can be shown (see [lB| 
for details) that the associated rate function iJ(/Lt) = 
limt_,oo(l/i) lnp(/i,t) satisfies 

H{^i) ~ DiSt-^'^h{Si^'^p./Di) for St > 1. (51) 

This scaling is confirmed numerically (inset of Fig. [T^ . 
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A ~ cDiSt^^/s St > 1 , (50) 

where c is a parameter-independent positive constant. 
Note that the original derivation [4^ of this law applies 
also to compressible carrier flows, so the constant c de- 
pends on the compressibility of the fluid velocity fleld. It 



FIG. 12: Lyapunov exponent A versus St. The dashed line is 
the asymptotic prediction (|50|l . Inset: rate function H{^) for 
various large values of St. 

We finally comment on how the stretching rate fluctu- 
ations change with St. Taylor expansion of H around its 
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minimum together with the scaling behavior (jSip shows 
that the standard deviation of the stretching rate is of 
the order of St^^^^ / \/t. For a given time i, the stretch- 
ing rate /i distributes more and more sharply around A 
when St increases. This behavior was anticipated by the 
numerical measurements reported in Section IIVI and is 
observed in direct numerical simulations of heavy parti- 
cles in homogeneous isotropic flows [27i] . 



VIII. REMARKS AND CONCLUSIONS 

Before concluding this paper the results discussed so 
far are commented in the light of what is known about 
real turbulent suspensions, which are relevant to most 
applications. Let us start by recalling the main fea- 
tures of turbulent flows. Turbulence is a multi-scale phe- 
nomenon [4^ which spans length scales ranging from a 
large (energy injection) scale L to the very small (dissi- 
pative) scale 77, often called the Kolmogorov scale. This 
hierarchy of length scales is associated with a hierarchy 
of time scales: from the large-scale eddy turnover time tl 
to the Kolmogorov time r^. Both ratios L/r] and tl/t,-i 
increase with the Reynolds number Re of the turbulent 
flow. Therefore, in general settings, no separation of time 
scales can be invoked to simplify the motion of suspended 
particles. However, in two circumstances simplifications 
are possible, namely: 

(i) For particles with a response time t much greater 
than T/^, the fluid velocity seen by the particle can be 
approximated by a random flow belonging to the Kraich- 
nan ensemble, as discussed in this paper. Then a Holder 
exponent /i = 1 or /i < 1 is chosen to study the dissipa- 
tive or inertial scales of turbulence, respectively. 

(ii) For intermediate response times ^ r ^ tl, at 
least for single or two-particle motions, the fluid veloc- 
ity seen by the particles can be approximated by an 
anisotropic generalization of the Kraichnan model }47| . 

In both asymptotics, the Kraichnan model and its 
generalization allow for predictions on single- and 
two-particle properties, many of them were discussed 
throughout this paper. In the following we discuss 
them in the context of turbulent suspensions. We fo- 
cus mostly on two-particle properties at dissipative and 
inertial scales. 

Dissipative range At such small scales, particles form 
(multi)fractal clusters, which can be quantitatively char- 
acterized by the St-dependence of the correlation dimen- 
sion I?2 or, equivalently, of the dimensional deficit d— 1?2 
(in turbulence one can define St = r/T^). Numerical stud- 
ies [l3, HO] show that the qualitative St-dependence of 
I?2 is similar to that observed in the Kraichnan model. 
Despite such similarities, it is likely that in turbulence, 
ejection from vortical regions play, at least for small St, 
an important role f50| . This can clearly not be accounted 
for in Kraichnan flows, as ^-correlated flows have no per- 
sistent structures. The absence of time correlations cer- 
tainly affects also the scaling behavior when St ^ 1 of the 



dimension deficit: while in turbulence [5|, |40| and time- 
correlated stochastic flows [ill . [39| it is observed that 
d— 1?2 Of St^, we have shown here that the behavior is 
linear in St. These discrepancies originate from the fact 
that white-in-time carrier flows are valid approximations 
of turbulence only for St3>l. 

Another question concerns the relative dispersion of 
a particle pair. In the dissipative range, the velocity 
field is smooth, so that particles separate exponentially 
with a rate given by the largest Lyapunov exponent A. 
If r ^ tl the results presented in previous Sections 
should apply, i.e. A oc St~^/^. For r,, <C r <C r^, the 
anisotropic generalization of the Kraichnan model pre- 
dicts A c!c St~^/^ [13] ■ However, the measurements of 
Lyapunov exponents made up to now (see e.g. [S^]) do 
not involve high-enough Stokes and Reynolds numbers to 
test the validity of these predictions in turbulent flows 

Inertial range As shown in this paper, for rough 
Kraichnan-type carrier flows, particles also form clusters 
which are however not fractal as they were in the dissi- 
pative range. This seems to be in qualitative agreement 
with the observations made in the inertial range of turbu- 
lence: Inhomogeneities have been found in 2d turbulence 
in the inverse cascade regime [i^, Eoj as well as in 3(i 
turbulence [sol . [5l[. However, while in the Kraichnan 
case the particle distribution depends on the local Stokes 
number St(r) only, this does not seem to be the case in 
turbulence, at least for St(r) <C 1 as studied in [HO] ( 
which in turbulence is defined by St(r) ^r/r^, being 
the characteristic turbulent time scale associated to the 
scale r). In turbulent flows, for small values of St(r), 
a different rescaling related to that of the acceleration 
(and hence pressure) field has been found [13]. How- 
ever such discrepancies do not question the relevance of 
the Kraichnan model to turbulent flows as it is expected 
to be a good approximation only for scales r such that 
Tr <^ r, i.e. St(r) 3> 1. Experiments or direct numerical 
simulations with high Re and St are thus needed to actu- 
ally test the validity of the dynamical scaling in terms of 
St(r) and to reproduce an equivalent of Fig. [3] for turbu- 
lent flows. As far as particle separation is concerned, we 
have seen in Section HVl that at very long times, and thus 
for separations r such that t Tr one should expect to 
observe Richardson dispersion. For intermediate times 
at which the sep aration is such that ^ ^ t, it is 
predicted in [47'] that an intermediate asymptotic regime 
may emerge with the typical particle separation r grow- 
ing as i^, i.e. much faster than Richardson diffusion. On 
the numerical and experimental side, we are not aware of 
any results on the relative dispersion of two heavy parti- 
cles in the inertial range. Testing the above predictions 
can be probably done only in experiments where Re can 
be very high. 

In summary, this paper reviews most of current under- 
standing of heavy particle suspensions in Kraichnan-like 
stochastic flows. In particular, we examined in details 
two-particle statistics both in smooth and rough veloc- 
ity fields. Numerical simulations, validated by analytics 
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originally derived in this paper, show that particle clus- 
tering is more efficient for smooth than rough flows, and 
can be characterized in terms of the local Stokes num- 
ber. Detailed predictions can be done in the very small 
and very large Stokes number asymptotics. In the former 
we provided an analytical expression for the dimensional 
deficit for any value of the fluid Holder exponent. More 
specifically, it is shown that the departure from a uni- 
form distribution is linear in the Stokes number, a result 
which is confirmed by numerics. As for the evolution 
of the relative separation of particle pairs at small sepa- 
rations, a well-verified asymptotic behavior for the Lya- 
punov exponent is discussed. At larger scales, by con- 
verting the scale-dependent Stokes number into a time- 
dependent one, we provided an original way to account 
for the recovering of tracer-like Richardson diffusion. Fi- 
nally, the relevance of these results, together with other 
predictions obtained in recent years from Kraichnan-like 
models of heavy particle suspensions, to particles in tur- 
bulent flows has been discussed. 

To conclude this work we suggest two different direc- 
tions for further investigations. First most of the pre- 



dictions related to the large-Stokes asymptotics lack nu- 
merical or experimental evidence in fluid flows with high 
Reynolds numbers and particles with huge inertia. Sec- 
ond it is now definitely clear that an important challenge 
for the near future is to understand whether or not some 
of the techniques developed for suspensions in random 
time-uncorrelated flows can be generalized/extended to 
time-correlated flows. For instance, a quantitative un- 
derstanding of the small-Stokes-number asymptotics in 
models that are closer to turbulence would be of great 
interest to many applications. A first step in this direc- 
tion has been recently attempted in [52 1. 
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